Let us set some global options for all code chunks in this document.
# Set seed for reproducibility
set.seed(1982)
# Set global options for all code chunks
knitr::opts_chunk$set(
# Disable messages printed by R code chunks
message = FALSE,
# Disable warnings printed by R code chunks
warning = FALSE,
# Show R code within code chunks in output
echo = TRUE,
# Include both R code and its results in output
include = TRUE,
# Evaluate R code chunks
eval = TRUE,
# Enable caching of R code chunks for faster rendering
cache = FALSE,
# Align figures in the center of the output
fig.align = "center",
# Enable retina display for high-resolution figures
retina = 2,
# Show errors in the output instead of stopping rendering
error = TRUE,
# Do not collapse code and output into a single block
collapse = FALSE
)# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(grateful)
library(ggplot2)
library(reshape2)
library(plotly)We want to solve the fractional diffusion equation \[\begin{equation} \label{eq:maineq} \partial_t u+(\kappa^2-\Delta_\Gamma)^{\frac{\alpha}{2}} u=f \text { on } \Gamma \times(0, T), \quad u(0)=u_0 \text { on } \Gamma, \end{equation}\] where \(u\) satisfies the Kirchhoff vertex conditions \[\begin{equation} \label{eq:Kcond} \left\{\phi\in C(\Gamma)\;\Big|\; \forall v\in V: \sum_{e\in\mathcal{E}_v}\partial_e \phi(v)=0 \right\} \end{equation}\]
If \(f=0\), then the solution is given by \[\begin{equation} \label{eq:sol_reprentation} u(s,t) = \displaystyle\sum_{j\in\mathbb{N}}e^{-\EIGENVAL^{\frac{\alpha}{2}}_jt}\left(u_0, e_j\right)_{L_2(\Gamma)}e_j(s). \end{equation}\]
# Function to build a tadpole graph and create a mesh
gets_graph_tadpole <- function(h){
edge1 <- rbind(c(0,0),c(1,0))
theta <- seq(from=-pi,to=pi,length.out = 10000)
edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
edges = list(edge1, edge2)
graph <- metric_graph$new(edges = edges)
graph$build_mesh(h = h)
return(graph)
}Let \(\Gamma_T = (\mathcal{V},\mathcal{E})\) characterize the tadpole graph with \(\mathcal{V}= \{v_1,v_2\}\) and \(\mathcal{E}= \{e_1,e_2\}\) as specified in Figure \(\ref{Interval.Circle.Tadpole}\)c. The left edge \(e_1\) has length 1 and the circular edge \(e_2\) has length 2. As discussed in Subsection \(\ref{subsec:prelim}\), a point on \(e_1\) is parameterized via \(s=\left(e_1, t\right)\) for \(t \in[0,1]\) and a point on \(e_2\) via \(s=\left(e_2, t\right)\) for \(t\in[0,2]\). One can verify that \(-\Delta_\Gamma\) has eigenvalues \(0,\left\{(i \pi / 2)^2\right\}_{i \in \mathbb{N}}\) and \(\left\{(i \pi / 2)^2\right\}_{2 i \in \mathbb{N}}\) with corresponding eigenfunctions \(\phi_0\), \(\left\{\phi_i\right\}_{i \in \mathbb{N}}\), and \(\left\{\psi_i\right\}_{2 i \in \mathbb{N}}\) given by \(\phi_0(s)=1 / \sqrt{3}\) and \[\begin{equation*} \phi_i(s)=C_{\phi, i}\begin{cases} -2 \sin (\frac{i\pi}{2}) \cos (\frac{i \pi t}{2}), & s \in e_1, \\ \sin (i \pi t / 2), & s \in e_2, \end{cases}, \quad \psi_i(s)=\frac{\sqrt{3}}{\sqrt{2}} \begin{cases} (-1)^{i / 2} \cos (\frac{i \pi t}{2}), & s \in e_1, \\ \cos (\frac{i \pi t}{2}), & s \in e_2, \end{cases}, \end{equation*}\] where \(C_{\phi, i}=1\) if \(i\) is even and \(C_{\phi, i}=1 / \sqrt{3}\) otherwise. Moreover, these functions form an orthonormal basis for \(L_2(\Gamma_T)\).
# Function to compute the eigenfunctions
tadpole.eig <- function(k,graph){
x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2])
x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2])
if(k==0){
f.e1 <- rep(1,length(x1))
f.e2 <- rep(1,length(x2))
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2)
f.e2 <- sin(pi*k*x2/2)
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
if((k %% 2)==1){
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
f.e2 <- cos(pi*k*x2/2)
f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1])
f <- list(phi=f1,psi=f2/sqrt(3/2))
}
}
return(f)
}Implementation of \(u\)
h <- 0.001
graph <- gets_graph_tadpole(h = h)
# Compute the FEM matrices
graph$compute_fem()
x <- graph$mesh$V[, 1]
y <- graph$mesh$V[, 2]
edge_number <- graph$mesh$VtE[, 1]
pos <- sum(edge_number == 1)+1
order_to_plot <- function(v)return(c(v[1], v[3:pos], v[2], v[(pos+1):length(v)], v[2]))
weights <- graph$mesh$weightskappa <- 1
n_finite <- 100
EIGENVAL <- c() # initialize empty vector for eigenvalues
EIGENFUN <- NULL # initialize NULL for eigenfunctions matrix
INDEX <- c()
for (j in 0:n_finite) {
lambda_j <- kappa^2 + (j*pi/2)^2
e_j <- tadpole.eig(j,graph)$phi#/sqrt(1000)
EIGENVAL <- c(EIGENVAL, lambda_j) # append scalar to vector
EIGENFUN <- cbind(EIGENFUN, e_j) # append column to matrix
INDEX <- c(INDEX, j)
if (j>0 && (j %% 2 == 0)) {
lambda_j <- kappa^2 + (j*pi/2)^2
e_j <- tadpole.eig(j,graph)$psi#/sqrt(1000)
EIGENVAL <- c(EIGENVAL, lambda_j) # append scalar to vector
EIGENFUN <- cbind(EIGENFUN, e_j) # append column to matrix
INDEX <- c(INDEX, j+0.1)
}
}
# Compute inner product matrix
INNER_PRODUCTS <- t(EIGENFUN) %*% (EIGENFUN * weights)
# Assuming INNER_PRODUCTS is a numeric matrix
plot_ly(
z = INNER_PRODUCTS,
type = "heatmap",
colorscale = "Viridis" # You can try "Hot", "Jet", etc.
) %>%
layout(
title = "Inner Product Matrix of Eigenfunctions",
xaxis = list(title = "Index"),
yaxis = list(title = "Index")
)EIGENFUN <- apply(EIGENFUN, 2, order_to_plot)
x <- order_to_plot(x)
y <- order_to_plot(y)
# Create interactive plot
fig <- plot_ly(x = ~INDEX, y = ~EIGENVAL, type = 'scatter', mode = 'markers',
marker = list(size = 6),
name = "Eigenvalues")
fig <- fig %>% layout(title = "Eigenvalues",
xaxis = list(title = "index i"),
yaxis = list(title = "Eigenvalues_i"))
plot_data <- data.frame(
x = rep(x, times = ncol(EIGENFUN)),
y = rep(y, times = ncol(EIGENFUN)),
z = as.vector(EIGENFUN),
frame = rep(INDEX, each = length(x))
)
# Compute axis limits
x_range <- range(x)
y_range <- range(y)
z_range <- range(EIGENFUN)
# Initial plot setup (first frame only)
p <- plot_ly(plot_data, frame = ~frame) %>%
add_trace(
x = ~x, y = ~y, z = ~z,
type = "scatter3d", mode = "lines",
name = paste0("n_finite = ", n_finite),
line = list(color = "blue", width = 2)
) %>%
layout(
scene = list(
xaxis = list(title = "x", range = x_range),
yaxis = list(title = "y", range = y_range),
zaxis = list(title = "Value", range = z_range)
),
updatemenus = list(
list(
type = "buttons", showactive = FALSE,
buttons = list(
list(label = "Play", method = "animate",
args = list(NULL, list(frame = list(duration = 100, redraw = TRUE), fromcurrent = TRUE))),
list(label = "Pause", method = "animate",
args = list(NULL, list(mode = "immediate", frame = list(duration = 0), redraw = FALSE)))
)
)
),
title = "Time: 0"
)
# Convert to plotly object with frame info
pb <- plotly_build(p)
# Inject custom titles into each frame
for (i in seq_along(pb$x$frames)) {
t <- INDEX[i]
err <- signif(EIGENVAL[i], 4)
pb$x$frames[[i]]$layout <- list(title = paste0("i: ", t, " | Eigenvalue_i: ", err))
}
fig # Display the plot